# R script for paper Manger - Pickup
# The Coevolution of Trade Agreement Networks and Democracy
# Journal of Conflict Resolution

#setwd("pathname /Replication Files")

# read data
net94 <- as.matrix(read.table("CoevPTA1994.dat"))
net95 <- as.matrix(read.table("CoevPTA1995.dat"))
net96 <- as.matrix(read.table("CoevPTA1996.dat"))
net97 <- as.matrix(read.table("CoevPTA1997.dat"))
net98 <- as.matrix(read.table("CoevPTA1998.dat"))
net99 <- as.matrix(read.table("CoevPTA1999.dat"))
net00 <- as.matrix(read.table("CoevPTA2000.dat"))
net01 <- as.matrix(read.table("CoevPTA2001.dat"))
net02 <- as.matrix(read.table("CoevPTA2002.dat"))
net03 <- as.matrix(read.table("CoevPTA2003.dat"))
net04 <- as.matrix(read.table("CoevPTA2004.dat"))

democ  <- as.matrix(read.table("udsInt94-04_1.dat"))
lowinc <- as.matrix(read.table("LowIncDum94-04.dat"))
medinc <- as.matrix(read.table("MidIncDum94-04.dat"))
hiinc  <- as.matrix(read.table("HiIncDum94-04.dat"))
invGDP <- as.matrix(read.table("invGDP94-04.dat"))

dist <- as.matrix(read.table("Coevdist.dat"))

Ally94  <- as.matrix(read.table("CovAlly1994.dat"))
Ally95  <- as.matrix(read.table("CovAlly1995.dat"))
Ally96  <- as.matrix(read.table("CovAlly1996.dat"))
Ally97  <- as.matrix(read.table("CovAlly1997.dat"))
Ally98  <- as.matrix(read.table("CovAlly1998.dat"))
Ally99  <- as.matrix(read.table("CovAlly1999.dat"))
Ally00  <- as.matrix(read.table("CovAlly2000.dat"))
Ally01  <- as.matrix(read.table("CovAlly2001.dat"))
Ally02  <- as.matrix(read.table("CovAlly2002.dat"))
Ally03  <- as.matrix(read.table("CovAlly2003.dat"))


Trade94  <- as.matrix(read.table("Cvtrade1994.dat"))
Trade95  <- as.matrix(read.table("Cvtrade1995.dat"))
Trade96  <- as.matrix(read.table("Cvtrade1996.dat"))
Trade97  <- as.matrix(read.table("Cvtrade1997.dat"))
Trade98  <- as.matrix(read.table("Cvtrade1998.dat"))
Trade99  <- as.matrix(read.table("Cvtrade1999.dat"))
Trade00  <- as.matrix(read.table("Cvtrade2000.dat"))
Trade01  <- as.matrix(read.table("Cvtrade2001.dat"))
Trade02  <- as.matrix(read.table("Cvtrade2002.dat"))
Trade03  <- as.matrix(read.table("Cvtrade2003.dat"))

# Prepare data for RSiena.
library(RSiena)

# create dependent variable
PTAs <- sienaNet(array(c(net94, net95, net96, net97, net98,
               net99, net00, net01, net02, net03, net04),
                        dim = c(145, 145, 11) ))

#create dependent behavior variable 
Democracy <- sienaNet(democ, type="behavior")

# create monadic time varying covariates
GDP_inv   <- varCovar(invGDP)

# create constant dyadic covariate
Distance  <- coDyadCovar(dist)

# create dyadic time varying covariates
Alliance  <- varDyadCovar(array(c(Ally94, Ally95, Ally96, Ally97, Ally98,
                                   Ally99, Ally00, Ally01, Ally02, Ally03),
                          dim = c(145, 145, 10) ))

Trade   <- varDyadCovar(array( c(Trade94, Trade95, Trade96, Trade97, Trade98,
               Trade99, Trade00, Trade01, Trade02, Trade03),
                        dim = c(145, 145, 10) ))

# define data                          
PTAbehdata <- sienaDataCreate(PTAs, Democracy, GDP_inv, Distance,
                           Alliance, Trade)
# create effects structure
PTAeff <- getEffects(PTAbehdata )

# give global report about data
print01Report(PTAbehdata, PTAeff, modelname = 'PTA_coev9404_UDS_int' )
   
# Define the kind of non-directed model and an extra long phase 3.
model33 <- sienaModelCreate(projname = "PTA_coev9404_UDS_int", modelType = 3, n3=9000)

# Specify the model
PTAeff <- setEffect(PTAeff, transTriads, include=FALSE)
PTAeff <- includeEffects(PTAeff, nbrDist2)
PTAeff <- includeEffects(PTAeff, X, interaction1="Distance")
PTAeff <- includeEffects(PTAeff, X, interaction1="Alliance")
PTAeff <- includeEffects(PTAeff, altX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="Democracy")
PTAeff <- includeEffects(PTAeff, altX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, egoXaltX, interaction1="GDP_inv")
PTAeff <- includeEffects(PTAeff, X, interaction1="Trade")
PTAeff <- includeInteraction(PTAeff, X, altX, interaction1=c("Trade", "GDP_inv"))
PTAeff <- includeEffects(PTAeff, name = "Democracy", avAlt, outdeg, interaction1 = "PTAs")
PTAeff <- setEffect(PTAeff, name = "Democracy", quad, include=TRUE)

# Check the current model specification
PTAeff

# estimate
ans3 <- siena07(model33, data = PTAbehdata, effects = PTAeff)
# alternatively if your computer permits, you can use parallel processing in R
# to speed up the simulations.
# library(parallel)
# ans3 <- siena07(model33, data = PTAbehdata, effects = PTAeff, batch=TRUE, nbrNodes = 4, useCluster = TRUE)
ans3


